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Abstract 

In this paper, we develop efficient randomized algorithms for estimating probabilistic robust- 
ness margin and constructing robustness degradation curve for uncertain dynamic systems. One 
remarkable feature of these algorithms is their universal applicability to robustness analysis prob- 
lems with arbitrary robustness requirements and uncertainty bounding set. In contrast to existing 
probabilistic methods, our approach does not depend on the feasibility of computing deterministic 
robustness margin. We have developed efficient methods such as probabilistic comparison, prob- 
abilistic bisection and backward iteration to facilitate the computation. In particular, confidence 
interval for binomial random variables has been frequently used in the estimation of probabilistic 
robustness margin and in the accuracy evaluation of estimating robustness degradation function. 
Motivated by the importance of fast computing of binomial confidence interval in the context 
of probabilistic robustness analysis, we have derived an explicit formula for constructing the 
confidence interval of binomial parameter with guaranteed coverage probability. The formula 
overcomes the limitation of normal approximation which is asymptotic in nature and thus in- 
evitably introduce unknown errors in applications. Moreover, the formula is extremely simple 
and very tight in comparison with classic Clopper-Pearson's approach. 

1 Introduction 

In recent years, there have been growing interest on the development of probabilistic methods for 
robustness analysis and design problems aimed at overcoming the computational complexity and 
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the issue of conservatism of deterministic worst case framework [21], [22], [191 HI I2H 121 HI [231 El 
El [lOl [261 123 EQ] . In the deterministic worst case framework, one is interested in computing 
the deterministic robustness margin which is defined as the maximal radius of uncertainty set such 
that the robustness requirement is guaranteed everywhere. However, it should be borne in mind 
that the uncertainty set may include worst cases which never happen in reality. Instead of seeking 
the worst case guarantee, it is sometimes "acceptable" that the robustness requirement is satisfied 
for most of the cases. It has been demonstrated that the proportion of systems guaranteeing the 
robustness requirement can be close to 1 even if the radii of uncertainty set are much larger than the 
deterministic robustness margin. The idea of sacrificing extreme cases of uncertainty has become a 
new paradigm of robust control. In particular, the concept of probabilistic robustness margin has been 
recently established [21 |2H El SI El 123] • Moreover, to provide more insight for the robustness of the 
uncertain system, the concept of robustness degradation function has been developed by a number 
of researchers [21 E]. For example, Barmish and Lagoa [5] have constructed a curve of robustness 
margin amplification versus risk in a probabilistic setting. In a similar spirit, Calafiore, Dabbene 
and Tempo [61 [7] have constructed a probability degradation function in the context of real and 
complex parametric uncertainty. Such a function describes quantitatively the relationship between 
the proportion of systems guaranteeing the robustness requirement and the radius of uncertainty 
set. Clearly, it can serve as a guide for control engineers in evaluating the robustness of a control 
system once a controller design is completed. It is important to note that the robustness degradation 
function and the probabilistic robustness margin can be estimated in a distribution free manner. This 
can be justified by the Truncation Theory established by Barmish, Lagoa and Tempo [2] and can 
also be illustrated by relaxing the deterministic worst-case paradigm. 

Existing techniques for constructing the robustness degradation curve relies on the feasibility of 
computing the deterministic robustness margin. Obtaining the probabilistic robustness margin is 
possible only when the robustness degradation curve is constructed. However, the computation of 
deterministic robustness margin is possible only when the robustness requirement P is simple and the 
uncertainty bounding set is of some special structure. For example, P is stability requirement and 
uncertainty bounding sets are spectral normed balls. In that case, deterministic analysis methods 
such as /i-theory can be applied to compute the deterministic robustness margin. The determinis- 
tic robustness margin is then taken as a starting point of uncertainty radius interval for which the 
robustness degradation curve is to be constructed. In general, the problem of computing the deter- 
ministic robustness margin is not tractable. Hence, to construct a robustness degradation curve of 
practical interest, the selection of uncertainty radius interval is itself a question. Clearly, the range of 
uncertainty radius for which robustness degradation curve is significantly below 1 is not of practical 
interest since only a small risk can be tolerated in reality. From application point of view, it is only 
needed to construct robustness degradation curve for the range of uncertainty radius such that the 
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curve is above an a priori specified level 1 — e where risk parameter e € (0, 1) is acceptably small. 

In this work, we consider robustness analysis problems with arbitrary robustness requirement and 
uncertainty bounding set. We develop efficient randomized algorithms for estimating probabilistic 
robustness margin p e which is defined as the maximal uncertainty radius such that the probability 
of guaranteeing the robust requirements is at least 1 — e. We have also developed fast algorithms for 
constructing robustness degradation curve which is above a priori specified level 1 — e. In particular, 
we have developed efficient mechanisms such as probabilistic comparison, probabilistic bisection and 
backward iteration to reduce the computational complexity. Complexity of probabilistic comparison 
techniques are rigorously quantified. 

In our algorithms, confidence interval for binomial random variables has been frequently used to 
improve the efficiency of estimating probabilistic robustness margin and in the accuracy evaluation 
of robustness degradation function. Obviously, fast construction of binomial confidence interval is 
important to the efficiency of the randomized algorithm. Therefore, we have derived an explicit 
formula for constructing the confidence interval of binomial parameter with guaranteed coverage 
probability [11] . The formula overcomes the limitation of normal approximation which is asymptotic 
in nature and thus inevitably introduce unknown errors in applications. Moreover, the formula is 
extremely simple and very tight in comparison with classic Clopper-Pearson's approach. 

The paper is organized as follows. Section 2 is the problem formulation. Section 2 discusses 
binomial confidence interval. Section 3 is devoted to probabilistic robustness margin. Section 4 
presents algorithms for constructing robustness degradation curve. Illustrative examples are given 
in Section 5. Section 6 is the conclusion. Proofs are included in Appendix. 

2 Problem Formulations 

We adopt the assumption, from the classical robust control framework, that the uncertainty is 
deterministic and bounded. We formulate a general robustness analysis problem in a similar way as 
[TO] as follows. 

Let P denote a robustness requirement. The definition of P can be a fairly complicated combi- 
nation of the following: 

• Stability or P-stability; 

• norm of the closed loop transfer function; 

• Time specifications such as overshoot, rise time, settling time and steady state error. 

Let B(r) denote the set of uncertainties with size smaller than r. In applications, we are usually 
dealing with uncertainty sets such as the following: 
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• Ip ball B p (r) := {A E R n : ||A|| P < r} where \\.\ p denotes the l p norm and p = 1, 2, ■ • ■ , oo. In 
particular, B°°(r) denotes a box. 

• Spectral norm ball B a {r) := {A G A : cr(A) < r} where cr(A) denotes the largest singular 
value of A. The class of allowable perturbations is 

A := {blockdiag[gi/ ri , • • • , q s I rs , Ai, ■ • ■ ,A C ]} (1) 

where qi G F, i = 1, • • • , s are scalar parameters with multiplicity • • ■ ,r s and Aj E 
F n ' xmi , i = 1, •■• ,c are possibly repeated full blocks. Here F is either the complex field 
C or the real field R. 

• Homogeneous star-shaped bounding set Bh(t) '■= {r(A — Ao) + Ao : A € Q} where Q C R n 
and Ao € Q (see [2] for a detailed illustration). 

Throughout this paper, B(r) refers to any type of uncertainty set described above. Define a function 
£(.) such that, for any X, 

£(X) := min{r : X G B(r)}, 
i.e., B(£(X)) includes X exactly in the boundary. By such definition, 

£(X) = min (r : - - A ° + A G Q 



£(X)=a(X), 

and 

£(X) = \\X\\ p 

in the context of homogeneous star-shaped bounding set, spectral norm ball and l p ball respectively. 

To allow the robustness analysis be performed in a distribution-free manner, we introduce the 
notion of proportion as follows. For any A G B{r) there is an associated system G(A). Define 
proportion 

vol({A G B(r) : The associated system G(A) guarantees P}) 
(r) := vol(£(r)) 

with 



vol(S) := / dq, 
Jqes 

where the notion of dq is illustrated as follows: 



(I) : If q = [x rs ] nxm is a real matrix in R nxm , then dq = Ylr=i YiT=i dx rs . 

(II) : If q = [x rs + jUrslnxm is a complex matrix in C nxm , then dq = n™=l Y\T=i{dx r sdy r 
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• (III): If q G A, i.e., q possess a block structure denned by (pQ), then dq = (111=1 (Ili=i 
where the notion of dqi and dAj is defined by (I) and (II). 

It follows that P(r) is a reasonable measure of the robustness of the system [6j[25]- In the worst 
case deterministic framework, we are only interested in knowing if P is guaranteed for every A. 
However, one should bear in mind that the uncertainty set in our model may include worst cases 
which never happen in reality. Thus, it would be "acceptable" in many applications if the robustness 
requirement P is satisfied for most of the cases. Hence, due to the inaccuracy of the model, we should 
also obtain the value of P(r) for uncertainty radius r which exceeds the deterministic robustness 
margin. 

Clearly, P(r) is deterministic in nature. However, we can resort to a probabilistic approach to 
evaluate P(r). To see this, one needs to observe that a random variable with uniform distribution 
over B(r), denoted by A", guarantees that 

vol(£>(r)J 

for any S, and thus 

P(r) = Pr{The associated system G(A. U ) guarantees P}. 

Define a Bernoulli random variable X such that X takes value 1 if the associated system G(A") 
guarantees P and takes value otherwise. Then estimating P(r) is equivalent to estimating binomial 
parameter Fx '■= Pr{X = 1} = P(r). It follows that a Monte Carlo method can be employed to 
estimate P(r) based on i.i.d. observations of X. 

Obviously, the robustness problem will be completely solved if we can efficiently estimate P(r) 
for all r G (0, oo). However, this is infeasible from computational perspective. In practice, only a 
small risk e can be tolerated by a system. Therefore, what is really important to know is the value of 
P(r) over the range of uncertain radius r for which P(r) is at least 1 — e where e G (0, 1) is referred 
as the risk parameter in this paper. When the deterministic robustness margin 

po := sup {r : G(A) guarantees P for any A G B(r)} 

can be computed, we can construct the robustness degradation curve as follows: (I) Estimate pro- 
portion for a sequence of discrete values of uncertainty radius which is started from r = po and 
gradually increased; (II) The construction process is continued until the proportion is below 1 — e. 
This is the conventional approach and it is feasible only when computing po is tractable. Unfor- 
tunately, in general there exists no effective technique for computing the deterministic robustness 
margin p$. In light of this situation, we establish a new approach which does not depend on the 
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feasibility of computing the deterministic robustness margin. Our strategy is to firstly estimate the 
probabilistic robustness margin 



and consequently construct the robust degradation curve in a backward direction (in which r is 
decreased) by choosing the estimate of p(e) as the starting uncertainty radius. 

To reduce computational burden, the estimation of probabilistic robustness margin relies on the 
frequent use of binomial confidence interval. The confidence interval is also served as a validation 
method for the accuracy of estimating robustness degradation function. Hence, it is desirable to 
quickly construct binomial confidence interval with guaranteed coverage probability. 

3 Binomial Confidence Intervals 

Clopper and Pearson [12] provided a rigorous approach for constructing binomial confidence interval. 
However, the computational complexity involved with this approach is very high. The standard 
technique is to use normal approximation which is not accurate for rare events. The coverage 
probability of the confidence interval derived from normal approximation can be significantly below 
the specified confidence level even for very large sample size. In the context of robustness analysis, 
we are dealing with rare events because the probability that the robustness requirement is violated 
is usually very small. We shall illustrate these standard methods as follows. 

3.1 Clopper-Pearson Confidence Limits 

Let the sample size ./V and confidence parameter 5 € (0, 1) be fixed. We refer an observation of X 
with value 1 as a successful trial. Let K denote the number of successful trials during the N i.i.d. 
sampling experiments. Let A; be a realization of K. The classic Clopper-Pearson lower confidence 
limit Lpf^^s and upper confidence limit f/jv,fc,<5 are given respectively by 



p(e) := sup{r : P(r) > 1 - e} 




if k = 



if k > 



and 




where p € (0, 1) is the solution of the following equation 




(2) 



and p € (0, 1) is the solution of the following equation 




k 



(3) 
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3.2 Normal Approximation 

It is easy to see that the equations ([2]) and ([3]) are hard to solve and thus the confidence limits are 
difficult to determine using Clopper-Pearson's approach. For large sample size, it is computationally 
intensive. To get around the difficulty, normal approximation has been widely used to develop simple 
approximate formulas (see, for example, [TS] and the references therein). Let $(.) denote the normal 
distribution function and Z§_ denote the critical value such that §>(Zs) = 1 — f . It follows from the 

2 2 

Central Limit Theorem that, for sufficiently large sample size N, the lower and upper confidence 

~ b I — (1— — ) ~ h I —(!——) 

limits can be estimated respectively as L « -ft — Zs_ y N N N and U « -ft + Zs y N N N . 

The critical problem with the normal approximation is that it is of asymptotic nature. It is not 
clear how large the sample size is sufficient for the approximation error to be negligible. Such an 
asymptotic approach is not good enough for studying the robustness of control systems. 

3.3 Explicit Formula 

It is desirable to have a simple formula which is rigorous and very tight for the confidence interval 
construction. Recently, we have derived the following simple formula for constructing the confidence 
limits [11] (The proof is provided in Appendix for purpose of completeness). 

Theorem 1 Let C(N, k, S) = * + f " V 1+gjV N> and U(N, k, 5) = *. + f " V 1+gJV " J 

TOift = 9 . T/ien Pr {£(N, K, 5) <F X < U(N, K, 5)} > 1 - <J. 

81n s 

Figures [1] and [2] show the confidence limits derived by different methods (curve A and B represent 
respectively the upper and lower confidence limits computed by Theorem [lj curve C and D represent 
respectively the upper and lower confidence limits calculated by Clopper-Pearson's method). It can 
be seen from these figures that our formula is very tight in comparison with the Clopper-Pearson's 
approach. Obviously, there is no comparison on the computational complexity. Our formula is 
simple enough for hand calculation. Simplicity of the confidence interval is especially important in 
the context of our robustness analysis problem since the confidence limits are repeatedly used for a 
large number of simulations. 



4 Estimating Probabilistic Robustness Margin 

In this section, we shall develop efficient randomized algorithms for constructing an estimate for p(e). 
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Figure 1: Confidence Interval (N = 1000, 5 = 10~ 2 ) 




Number of Successful Trials 

Figure 2: Confidence Interval (N = 1000, 6 = 0.001) 
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4.1 Separable Assumption 

We assume that the robustness degradation curve of the system can be separated into two parts by a 
horizontal line with height 1 — e, i.e., 

P(r) < 1 - e for all r > p(e). 

We refer such an assumption as the Separable Assumption. From an application point of view, 
this assumption is rather benign. Our extensive simulation experience indicated that, for small risk 
parameter e, most control systems guarantee the separable assumption. It should be noted that it is 
even much weaker than assuming that P(r) is non-decreasing (See illustrative Figure [3]). Moreover, 
the non-increasing assumption is rather mild. This can be explained by a heuristic argument as 
follows. Let 

B F (r) := {A € B(r) : The associated system G(A) guarantees P}. 



Then 



and 



cflP(r) 1 



dr vol(S(r)) 



_ vol(g p (r)) 
1 ' vol(i3(r)) 

dvol{B F {r)} dvol{B(r)} 



P( 



r 



dr dr 



(4) 



Moreover, due to the constraint of robust requirement P, it is true that B (r+dr) \ B (r) is a subset 
of B(r + dr) \ B(r) and it follows that vol{£>(r)} increases (as r increases) faster than vol{£> p (r)}, 
i.e., dv °^d r < d vol ^f^ r ^ ■ Hence inequality dvo1 ^, < P(r) - vol ^^ r ^ is not hard to guarantee 
in the range of r such that P(r) is close to 1. It follows from equation @ that —^p- < can be 
readily satisfied in the case of small risk parameter e. 

It is interesting to note that our randomized algorithms for estimating p{e) presented in the sequel 
rely only on the following assumption: 

P(r) < 1 - e Vr G (p(e), 2 K ] \J{¥ : k < i < 0} (5) 

where k = [log 2 p(e)~\ ■ It can be seen that condition ([S]) is even weaker than the separable assumption. 

When condition ([5]) is guaranteed, an interval which includes p(e) can be readily found by starting 
from uncertainty radius r = 1 and then successively doubling r or cutting r in half based on the 
comparison of P(r) with 1 — e. Moreover, bisection method can be employed to refine the estimate 
for p(e). Of course, the success of such methods depends on the reliable and efficient comparison 
of P(r) with 1 — e based on Monte Carlo method. In the following subsection, we illustrate a fast 
method of comparison. 
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Figure 3: Illustration of Separable Assumption. The robustness degradation curve can be separated 
as the upper segment and lower segment by the dashed horizontal line with height 1 — e. In this 
example, the separable assumption is satisfied, while the non-increasing assumption is violated. 

4.2 Probabilistic Comparison 

In general, Fx can only be estimated by a Monte Carlo method. The conventional method is to 
directly compare |r with 1 — e where K is the number of successful trials during N i.i.d. sampling 
experiments. There are three problems with the conventional method. First, the comparison of ^ 
with 1 — e can be very misleading. Second, the sample size N is required to be very large to obtain 
a reliable comparison. Third, we don't know how reliable the comparison is. In this subsection, 
we present a new approach which allows for a reliable comparison with many fewer samples. The 
key idea is to compare binomial confidence limits with the fixed probability 1 — e and hence reliable 
judgement can be made in advance. 

Function name: Probabilistic-Comparison. 

Input: Risk parameter e and confidence parameter 5. 
Output: d = Probabilistic-Comparison (5,e). 
Step 1. Let d <— 0. 

Step 2. While d = do the following: 
• Sample X. 
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Figure 4: Complexity of Probabilistic Comparison. The horizontal axis represents 1 — Px- The 
vertical axis represents sample size. The solid line and the dash-dot line respectively show the 
average sample size and the 95%-quantile of the sample size. 

• Update N and K. 

• Compute lower confidence limit L and and upper confidence limit U by Theorem 1. 

• If U < 1 - e then let d < 1. If L > 1 - e then let d <- 1. 

The confidence parameter 5 is used to control the reliability of the comparison. A typical value 
is 5 = 0.01. The implication of output is as follows: d = 1 indicates that Px > 1 — e is true with 
high confidence; d = — 1 indicates that Px < 1 — e is true with high confidence. 

Obviously, the sample size is random in nature. For e = 5 = 0.01, we simulated the Probabilistic 
Comparison Algorithm identically and independently 100 times for different values of Px- We observe 
that, for each value of Px, the Probabilistic Comparison Algorithm makes correct judgement among 
all 100 simulations. Figure U] shows the average sample size and the 95%-quantile of the sample size 
estimated from our simulation. It can be seen from the figure that, as long as Px is not very close to 
1 — e, the Probabilistic Comparison Algorithm can make a reliable comparison with a small sample 
size. 

4.3 Computing Initial Interval 

Under the separable assumption, an interval [a, b] which includes p{e) can be quickly determined by 
the following algorithm. 



11 



Function name: Initial. 

Input: Risk parameter e and confidence parameter 8. 
Output: [a, b] = Initial e). 

Step 1. Let r «— 1. Apply Probabilistic-Comparison algorithm to compare P(l) with 1 — e. Let the 
outcome be d±. 

Step 2. If di = 1 then let d <— c?i and do the following: 

• While (f = 1 do the following: 

— Let r <— 2r. Apply Probabilistic-Comparison algorithm to compare P(r) with 1 — e. 
Let the outcome be d. 

• Let a <— | and 6 <— r. 

Step 3. If e?i = — 1 then let d <— d\ and do the following: 

• While d = — 1 do the following: 

— Let r <— |. Apply Probabilistic-Comparison algorithm to compare P(r) with 1 — e. 
Let the outcome be d. 

• Let a <— r and 6 <— 2r. 
4.4 Probabilistic Bisection 

Once an initial interval [a, b] is obtained, an estimate R for the probabilistic robustness margin p(e) 
can be efficiently computed as follows. 

Function name: Bisection. 

Input: Risk parameter e, confidence parameter 5, initial interval [a, b], relative tolerance 7. 

Output: [R] = Bisection(7, a, b, 5, e). 

Step 1. Input 7, a, b, 5, e. 

Step 2. While b — a > 7a do the following: 

• Let r <— ^T^. Apply Probabilistic-Comparison algorithm to compare P(r) with 1 — e. Let 
the outcome be d. 

• If d = — 1 then let 6 <— r, else let a <— r. 
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Step 3. Return R = b (Note: this is actually a soft upper bound). 

It should be noted that when applying bisection algorithm to refine the initial interval [a, b], the 
execution of the algorithm may take very long time if P(r) ~ e. However, such chance is almost 0. 
This problem can be fixed by the following methods. 

• We can limit the maximum number of simulations to a number M. When conducting simulation 
at radius r, simulation results can be saved for uncertainty radius where a is the lower 
bound of the current interval with middle point r (See the next section for the idea of sample 
reuse). After the number of simulations for uncertainty radius r exceeds M, the simulation is 
switched for uncertainty radius ^f-- 

• In application, one might want to construct the robustness degradation curve in the backward 
direction by starting from an upper bound of p(e). In that situation, we don't need to compute 
a very tight interval for p(e) and hence the chance of P(r) ~ e is even smaller. 

5 Constructing Robustness Degradation Curve 

We shall develop efficient randomized algorithms for constructing robustness degradation curve, 
which provide more insight for the robustness of the uncertain system than probabilistic robustness 
margin. First we recall the Sample Reuse algorithm [10] for constructing robustness degradation 
curve for a given range of uncertainty radius. 

Sample Reuse Algorithm 

Input: Sample size N, confidence parameter 5 € (0, 1), uncertainty radius interval [a, 6], number of 
uncertainty radii I. 

Output: Proportion estimate Pj and the related confidence interval for r j = b — ( fe 1 ) ; j = 
1, 2, • • • , I. In the following, run denotes the number of sampling experiments conducted at r% 
and 7TT.j2 denotes the number of observations guaranteeing P during the mn sampling experi- 
ments. 

Step 1 (Initialization). Let M = [77i«]; X 2 be a zero matrix. 
Step 2 (Backward Iteration). For i = 1 to % = I do the following: 

• Let r <— Ti. 

• While run < N do the following: 

— Generate uniform sample q from B(r). Evaluate the robustness requirement P for q. 
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— Let m s i <— m s i + 1 for any s such that r s > £(q). 

— If robustness requirement P is satisfied for q then let m S 2 <— m S 2 + 1 for any s such 
that r s > t{q). 

• Let Pj <— and construct the confidence interval of confidence level 100(1 — 8)% by 
Theorem 1. 



It should be noted that the idea of the Sample Reuse Algorithm is not simply a save of sample 
generation. It is actually a backward iterative mechanism. In the algorithm, the most important 
save of computation is usually the evaluation of the complex robustness requirements P (See [10] for 
details). 

Now we introduce the global strategy for constructing robustness degradation curve. The idea is 
to successively apply the Sample Reuse Algorithm for a sequence of intervals of uncertainty radius. 
Each time the size of interval is reduced by half. The lower bound of the current interval is defined to 
be the upper bound of the next consecutive interval. The algorithm is terminated once the robustness 
requirement P is guaranteed for all N samples of an uncertainty set of which the radius is taken as 
the lower bound of an interval of uncertainty radius. More precisely, the procedure is presented as 
follows. 

Global Strategy 

Input: Sample size N, risk parameter e and confidence parameter 8 £ (0, 1). 
Output: Proportion estimate Pj and the related confidence interval. 
Step 1. Compute an estimate R for probabilistic robustness margin p(e). 
Step 2. Let STOP <- 0. Let a <- f and b <- R. 
Step 3 (Backward Iteration). While STOP = do the following: 

• Apply Sample Reuse Algorithm to construct robustness degradation curve for uncertainty 
radius interval [a,b]. 

• If the robustness property P is guaranteed for all iV samples of uncertainty set B(r) with 
radius r = a then let STOP <— 1, otherwise let b <— a and a *— k. 



For given risk parameter e and confidence parameter 5, the sample size is chosen as 



2(l-e+f )(1- f)ln| 
N > 6 
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with a E (0, 1). It follows from Massart's inequality [18] that such a sample size ensures 

K 



Pr 



Fx N 



< ae > > 1 



with Fx = 1 — e (See also Lemma Q] in Appendix). It should be noted that Massart's inequality is 
less conservative than the Chernoff bounds in both multiplicative and additive forms. 

We would like to remark that the algorithms we propose for estimating the probabilistic robust- 
ness margin and constructing robustness degradation curve are susceptible of further improvement. 
For example, the idea of sample reuse is not employed in computing the initial interval and in the 
probabilistic bisection algorithm. Moreover, in constructing the robustness degradation curve, the 
Sample Reuse Algorithm is independently applied for each interval of uncertainty radius. Actually, 
the simulation results can be saved for the successive intervals. 

We would also like to note that, for a practitioner, computing the probabilistic robustness margin 
might be sufficient for understanding the system robustness. However, when more insight about the 
system robustness is expected, the techniques introduced in Section 5 can be employed. Of course, 
the price is more computational effort. 



6 Illustrative Examples 

In this section we demonstrate through examples the power of our randomized algorithms in solving 
a wide variety of complicated robustness analysis problems which are not tractable in the classical 
deterministic framework. 

We consider a system which has been studied in [T3] by a deterministic approach. The system is 
as shown in Figure [5j 



C(s) 




P(s) 


c 


— ► 







Figure 5: Uncertain System 

The compensator is C(s) = j^jq and the plant is P(s) = s ( s+ 4^25^^+f>+o 3<5 3 ) w ^h parametric 
uncertainty A = [5±, 62, ^3] T - The nominal system is stable. The closeddoop roots of the nominal 
system are: 

Zl = -15.9178, z 2 = -1.8309, z 3 = -1.1256 + 7.3234i, z 4 = -1.1256 - 7.3234i. 
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The peak value, rise time, settling time of step response of the nominal system, are respectively, 
P peak = 1A7 > «? = 0.185, t° 8 = 3.175. 

We first consider the robust Instability of the system. The robustness requirement P is defined 
as P-stability with the domain of poles defined as: Real part < —1.5, or fall within one of the two 
disks centered at z% and with radius 0.3. The uncertainty set is defined as the polytope 

B H {r) := |rA + (l-r) ^ : A G conv{A 1 , A 2 , A 3 , A 4 } j 

where 'conv' denotes the convex hull of A* = [~ sin(^jp"7r), \ cos(^p 1 -7r), — -^] T for z = 1,2,3 and 

A 4 = [0, 0, 1] T 

Obviously, there exists no effective method for computing the deterministic robustness margin in 
the literature. However, our randomized algorithms can efficiently construct the robustness degra- 
dation curve. See Figure [6l 

In this example, the risk parameter is a priori specified as e = 0.001. The procedure for estimating 

the probabilistic robustness margin is explained as follows. Let N denote the sample size which is 

random in nature. Let K denote the number of successful trials among N i.i.d. sampling experiments 

as defined in Section 2 and Subsection 3.1 (i.e., a successful trial is equivalent to an observation that 

the robustness requirement is guaranteed). Let confidence parameter 5 = 0.01 and choose tolerance 

7 = 0.05. Staring from r = 1, after N = 7060 simulations we obtain K = 7060, the probabilistic 

comparison algorithm determined that P(l) > 1 — e since the lower confidence limit L > 1 — e. The 

simulation is thus switched to uncertainty radius r = 2. After N = 65 times of simulation, it is 

found that K = 61. The probabilistic comparison algorithm detected that P(2) < 1 — e because the 

upper confidence limit U < 1 — e. So, initial interval [1, 2] is readily obtained. Now the probabilistic 

bisection algorithm is invoked. Staring with the middle point of the initial interval (i.e., r = = |), 

after TV = 613 times of simulations, it is found that K = 607, the probabilistic comparison algorithm 

concluded that P(|) < 1 - e since the upper confidence limit U < 1 — e. Thus simulation is moved 

to r = — rp- = f • It is found that K = 9330 among N = 9331 times of simulations. Hence, the 

probabilistic comparison algorithm determined that P(|) > 1 — e since the lower confidence limit 

-+- 11 

L > 1 — e. Now the simulation is performed at r = = After iV = 6653 simulations, it is 
discovered that K = 6636. The probabilistic comparison algorithm judged that P(^) < 1 — e based 
on calculation that the upper confidence limit U < 1 — e. At this point the interval is [|, -^] and the 
bisection is terminated since tolerance condition b — a < 7a is satisfied. The evolution of intervals 
produced by the probabilistic bisection algorithm is as follows: 



[1.2] 







"5 3" 
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Now we have obtained an interval [|, which includes p(0.001), so the Sample Reuse Algorithm 
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can be employed to construct robustness degradation curve. In this example, the number of uncer- 
tainty radii is I = 100 and the confidence parameter is chosen as 5 = 0.001. A constant sample size 
is computed by formula (J6j) with a = 0.5 as 

N = 50,631. 

The interval from which we start constructing robustness degradation curve is ^-]. It is deter- 
mined that K = N = 50, 632 at uncertainty radius r = jq. Therefore, the Sample Reuse Algorithm 
is invoked only once and the overall algorithm is terminated (If K ^ N for r = j^, then the next 
interval would be ^]). Although P(r) is evaluated for I = 100 uncertainty radii with the same 
sample size N, the total number of simulations is only 153,358 << Nl = lOOiV. To provide an 
evaluation of accuracy for all estimates of P(r), confidence limits are computed by Theorem 1. 
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Figure 6: Robustness Degradation Curve 

We now apply our algorithms to a robustness problem with time specifications. Specifically, 
the robustness requirement P is : Stability, and rise time t r < 135% t® = 0.25, settling time 
t a < 110% t° s = 3.5, overshoot P peak < 116% P° eofe = 1.7. The uncertainty set is B^r) := {A : 
IIAHoo < r}. 

In this case, the risk parameter is a priori specified as e = 0.01. It is well known that, for this type 
of problem, there exists no effective method for computing the deterministic robustness margin in the 
literature. However, our randomized algorithms can efficiently construct the robustness degradation 
curve. See Figure [71 

We choose 7 = 0.25 and 5 = 0.01 for estimating p(0.01). Starting from uncertainty radius r = 1, 



Upper Confidence Limit 



Lower Confidence Limit 



Minimum Variance Unbias Estimate 
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the initial interval is easily found as [|, ^jr] through the following interval evolution: 



"1 

r\ 




"1 1" 




"i r 
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The sequence of intervals generated by the probabilistic bisection algorithm is as follows: 



"1 1" 




"1 3 " 




"1 5 " 


8' 4_ 




.8' 16. 
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So, we obtained an estimate for the probabilistic robustness margin p(0.01) as J^. To construct 
robustness degradation curve, the number of uncertainty radii is chosen as / = 100 and the confidence 
parameter is chosen as 5 = 0.01. A constant sample size is computed by formula © with a = 0.2 as 

N = 24,495. 

The interval from which we start constructing robustness degradation curve is [(g) Jj]- We found 
that this is also the last interval of uncertainty radius because it is determined that K = N at 
uncertainty radius r = Jr. 
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Figure 7: Robustness Degradation Curve 

7 Conclusions 

In this paper, we have established efficient techniques which applies to robustness analysis problems 
with arbitrary robustness requirements and uncertainty bounding set. The key mechanisms are 
probabilistic comparison, probabilistic bisection and backward iteration. Motivated by the crucial 
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role of binomial confidence interval in reducing the computational complexity, we have derived an 
explicit formula for computing binomial confidence limits. This formula overcomes the computational 
issue and inaccuracy of standard methods. 



A Proof of Theorem 1 

To show Theorem [JJ we need some preliminary results. The following Lemma [TJ is due to Massart 
US- 
Lemma 1 Pr {§ > F x + e} < exp (- 2(Px+| j V ( f_ Px _ |) ) for all e € (0, 1 - P x ). 

Of course, the above upper bound holds trivially for e > 1 — Fx- Thus, Lemma [T] is actually true 
for any e > 0. 

Lemma 2 Pr {§ < F x - e} < exp (- 2(J?x _^_ Px+ ^ for all e > 0. 

Proof. Define Y = 1 — X . Then Py = 1 — Fx- At the same time when we are conducting 
N i.i.d. experiments for X, we are also conducting N i.i.d. experiments for Y. Let the number 
of successful trials of the experiments for Y be denoted as Ky. Obviously, Ky = N — K . Ap- 
plying Lemma Q] to Y, we have Pr j^- > Fy + e j < exp y — 2 (p Y+ ±)^(i-t> Y -±) J ■ It follows that 

Pr { N ^j K > 1 — Fx + e} < exp ^— 2 (±-p x +c) [ i-(i-p x )-L ] J ■ The proof is thus completed by observ- 
ing that Pr > 1 - F x + e} = Pr {§ < F x - e}. □ 

The following lemma can be found in [13J . 
Lemma 3 X^=o i^)^ (^~ X ) N ~ : * decreases monotonically with respect to x € (0, 1) for k = 0, 1, • • • ,N. 

Lemma 4 £ - =0 {^(l-x) N ^ < exp (- - ( , x+ ^]~fK x _^ € (£, 1) for k = 0, 1, ■ ■ ■ ,N. 

Proof. Consider binomial random variable X with parameter Px > 77- Let K be the number of 
successful trials during N i.i.d. sampling experiments. Then X^=o (j )^x(I — F x ) N ~i = Pr{i"T < k}. 
Note that Pr{K < k} = Pr {§ < F x - {F x - 77)}- Applying Lemma [2] with e = F x - 77 > 0, we 
have 

-'*)-' * Jf'-*)' Fi _ A ) 

i=o^/ V 2(Px- ! ^)(l-Px + ? V^)/ 

= exp f N <F* ~ #) 2 ^ 
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Since the argument holds for arbitrary binomial random variable X with Fx > 77 , thus the proof of 
the lemma is completed. □ 

LemmaS E.to (^)^(l-x)^ > 1-exp f- ^-J^^-^ j Vx € (0, 77) for k = !,■■■ ,N. 

Proof. Consider binomial random variable X with parameter Fx < 77- Let K be the number of 
successful trials during N i.i.d. sampling experiments. Then 

E (Tj^xO- - ^x) N ~ J = <k} = Frl^<Fx + (^-F x )\. 

Applying Lemma Q] with e = 77 — Fx > 0, we have that 



}T V^ " V 2(P X + SZf2L) (1 _ Px _ 



fc m fc 



1 — exp 



iV(P x - |)2 



2 (§Px + &)(!- - sir: 



Since the argument holds for arbitrary binomial random variable X with Fx < 77 , thus the proof of 
the lemma is completed. □ 

Lemma 6 Let < k < N. Then L/v,fc,<5 < U^,k,s- 

Proof. Obviously, the lemma is true for k = 0, N. We consider the case that 1 < k < N — 1. 
Define S(N,k,x) := £*L 0^(1-^)^' for x G (0,1). Notice that S(N,k,p) = S(N,k- l,p) + 
(fc)p fc (l-p) W " fc = |- Thus 



s(#,fc-i,p)-s(j\r,ft-i,p) = i-^ 

Notice that 5 € (0, 1) and that p e (0, 1), we have that 



5(JV, fc - 1, p) - 5(iV, fc - 1, p) = 1 - S + ( ^ V(l - p) 7V ~ fc > 0. 



k 

By Lemma «S(iV, A; — 1, x) decreases monotonically with respect to x, we have that p < p. □ 
We are now in the position to prove Theorem [TJ For national simplicity, let 

P = U N , k) s, q = U(N,k,S). 

It can be easily verified that p < q for k = 0, iV. We need to show that p < q for < k < N. 
Straightforward computation shows that q is the only root of equation 



exp 



N(x - t>) 2 \_8 
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with respect to x € (-4, oo). There are two cases: q > 1 and q < 1. If g > 1 then p < g is trivially 
true. We only need to consider the case that < q < 1. In this case, it follows from Lemma E] that 

o— n \J / 



3=0 



Recall that 



we have 



2(|9+A)(l-§9-AV 2 ' 



3=0 

fc /ATN fc 



2' 



3=0 v J ' 3=0 

Therefore, by Lemma[3j we have that p < q for < k < N. Thus, we have shown that U^kS — ^N,k,s 
for all k. 

Similarly, by Lemma [5] and Lemma [3] we can show that g > C(N,k,5). By Lemma [BJ we 
have C(N,k,5) < Ljy t k,5 < UN,k,6 < U(N,k,S). Finally, the proof of Theorem CD is completed by 
invoking the probabilistic implication of the Clopper-Pearson confidence interval. 
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